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ABSTRACT 

We present a new definition of subhalos in dissipationless dark matter N-body simulations, based 
on the coherent identification of their dynamically bound constituents. Whereas previous methods 
of determining the energetically bound components of a subhalo ignored the contribution of all the 
remaining particles in the halo (those not geometrically or dynamically associated with the subhalo), 
our method allows for all the forces, both internal and external, exerted on the subhalo. We demon- 
strate, using the output of a simulation at different timesteps, that our new method is more accurate 
at identifying the bound mass of a subhalo. We then compare our new method to previously adopted 
means of identifying subhalos by applying each to a sample of 1838 virialized halos extracted from 
a high resolution cosmological simulation. We find that the subhalo distributions are similar in each 
case, and that the increase in the binding energy of a subhalo from including all the particles located 
within it is almost entirely balanced by the losses due to the external forces; the net increase in the 
mass fraction of subhalos is roughly 10%, and the extra substructures tending to reside in the inner 
parts of the system. Finally, we compare the subhalo populations of halos to the sub-subhalo popu- 
lations of subhalos, finding the two distributions to be similar. This is a new and interesting result, 
suggesting a self-similarity in the hierarchy substructures within cluster mass halos. 
Subject headings: cosmology: dark matter — galaxies: clusters: general — methods: N-body simula- 
tions 



1. INTRODUCTION 

In the current paradigm of hierarchical structure for- 
mation, satellite galaxies in clusters are associated with 
the remnants of dark matter halos - known as subhalos 
- that have, at some point in their history, been accreted 
and absorbed by a more massive halo. Once captured by 
their host (or mother) halo, subhalos are continuously 
eroded by the combined effects of dynamical friction and 
tidal stripping by the host halo core. Often a subhalo will 
lose a large fraction of its mass at the time of accretion, 
with only the dense core surviving until the present day. 
In studies of structure formation, it is these self-bound 
remnants of the original subhalo which we associate with 
satellite galaxies. 

With the advent of high resolution cosmological 
simulations, we have been able to test models of hi- 
erarchical structure formation in a ACDM universe. 
Until the end of the last decade, it was not possible 
to reach the required mass resolution to enable the 
identification of galaxy mass halos as substructures in 



cluste rs (IWhitdll976t Ivan Kampen| [l995: Summe rs et al.l 
Il995t iMoore et alj Il996f ). Commonly known as the 
overmerging problem, this was mainly due to the 
limited mass and force resolution of the simulations 
used. The major causes of this problem were premature 
tidal disruption due to inadequate force resolution and 
two-particle evaporation for ha los with a small number 
of particles (jKlypin et al.l Il999h . Over the last 7 years, 
rapid advances in parallel computing, through both 
the improvement of hardware and the development of 
fast and efficient parallel algorithms, has enabled us to 
achieve the numerical resolution required to overcome 
these numerical proble ms and probe the subhalo popula- 
tions of ACDM halo s dGhigna et all 119981: iKlvpin et al l 
[1991 IMoore et al l 119991: lOkamoto fc Habd 119991: 
Ghigna et all 120001: iBode et alj 120011; ISpringel et all 
boOlt iKravtsov "et al.l I2004U: iDe Lucia et alJ 120041 : 



iGao et al. 



2004k iReed et al.ll2005T) . 
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However, the exact definition of a 'subhalo' is some- 
what ambiguous, and it tends to be intrinsically linked 
to the algorithm adopted to identify the groups of par- 
ticles belonging to each subhalo in a numerical simu- 
lation. Therefore, what constitutes a substructure is 
often dependent on the process or quantity being in- 
vestigated. For example, if one is interested in using 
numerical simulations to investigate the impact of sub- 
structures within gravitational lenses, and therefore sam- 
pling the projected potential field within objects, then 
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a geometric al definition of a subhalo is sufficient (see , 
for examp le. Hcnnawi et al. (2005); Hag an et al.l (|2005[) ; 
iMao et all (|2004D ; IAmara et alj (12006ft ). However, if one 
instead wants to study the formation and evolution of 
subhalos as they initially grow through accretion, are 
captured by larger structures, and are subsequently de- 
pleted through dynamical friction and tidal ablation, one 
must then search for structures in 6 dimensional phase 
space in orde r to identify statistic a lly significant groups 
of particles dTaffoni et all 120031 lHayashi et al.l 120031: 
Kazantzid is et al.l l2004: Kra vtsov et a l. 2004b; Gill et all 



2004b). 



However, there is an additional and often dominant 
goal in defining 'substructure'. In most assemblages 
that we call 'galaxies', the lifetimes of the stars that 
we observe (especially in the K-band) are quite long 
compared to their dynamical (or orbital) times in the 
systems in which they are found. Hence, if the pur- 
pose is to identify objects that correspond to the galax- 
ies that we observe, then we must identify those sub- 
halo particles that will remain together over many dy- 
namical timescales. In this case, all the forces on 
a subhalo, internal and external, must be accounted 
for in order to identify such objects. Requiring that 
all particles be ene rgetically bound, as is adopted by 
many a uthors (e.g. iGhigna et all (12000ft: ISpringel et al.1 
(120011): iKravtsov et alJ (l2004bft: iDe Lucia et al l (12004ft: 
iGao et all (|2004D : iReed et all (|2005ft ). can be a neces- 
sary, but is not a sufficient, condition for coherence. We 
will later spell out what we believe to be a more appro- 
priate algorithm. 

There are many existing algorithms for defining and 
identifying halos and subhalos in dissipationless N-body 
simulations. As noted above, some methods are essen- 
tially geometrical, others aim at finding dynamically co- 
herent structures; all contain both dimcnsionless and di- 
mensional parameters. In general, there are three lev- 
els of refinement that one can adopt. The most basic 
is to use a geometrical rout i ne, such as the Friends- 
Of-Friends (iHuchra fc Geileil Il98i iDavis et all fl985l: 
Lacey fc Cole 1994ft or Denmax dBertschinger fc Gelbl 



19911: iGelb fe Bertschingerlll9 94: Eise nstein feHutll 1998ft 



algorithms. These use only instantaneous particle posi- 
tions to group together nearby particles (defined by the 
FOF 'linking length' or the Denmax 'smoothing length') 
into localized structures. Neither method performs any 
type of dynamical analysis to check whether these struc- 
tures are bound. 

The next level of refinement is to follow a geometrical 
identification of halos with a procedure for determining 
those particles that are dynamically associated with the 
substructure. Essentially, we wish to separate particles 
belonging to substructure from the local underlying back- 
ground matter belonging solely to the cluster. This we 
do in two stages, first removing the velocity outliers from 
the subhalo, and then refining the process by iteratively 
removing the unbound particles with the greatest total 
energy until only bound particles remain. This must be 
done iteratively to ensure that changes in the subhalo 
center of mass and velocity are accounted for as parti- 
cles are removed. Examples of publicly a vailable rou- 
tines that use this approach include SKID (IStadel et al.l 
fl997l iStadell l200l[ ). BDM (iKl vpin et al.l 11999ft. SUB- 
FIND (|Springel et all 12001ft . VOBOZ (iNevrinck et all 



[20051 ) and MHF (jGill et alj [2004a! ) . In each of these, 
the unbinding is performed assuming the subhalo is com- 
pletely isolated, i.e. only the bound particles in the sub- 
halo at each iteration are taken into consideration when 
calculating the potential energy. Not taken into account 
in this energy calculation are the (previously identified) 
unbound particles located spatially within the subhalo, 
nor the disruptive effect of the tidal forces from the par- 
ticles surrounding it. Therefore, in order to correctly 
identify halos and subhalos that are not only instanta- 
neously bound, but will also remain so in the near future 
(i.e. within the local dynamical timescale), one should 
account for all of the forces that may influence the dy- 
namical state of a subhalo. 

An instructive analogy is provided by the stars within 
the Milky Way. There are a significant number of stars 
that, by chance, have an extremely low velocity relative 
to that of the Sun. Consider one such star and the Sun as 
a single and isolated system, ignoring the effect of all the 
other matter in the Milky Way: by calculating the total 
energy of each star, it would appear that they are grav- 
itationally bound and therefore make up a binary pair. 
It is only once we also take into account all the other 
stars in the Milky Way that we realize that the gravi- 
tational force between the star in question and the Sun 
is dwarfed by the tidal forces due to all the other stars 
in the galaxy. It then becomes clear that the two stars 
do not constitute a stable and bound system. The same 
argument may be applied to dark matter substructures 
in a cluster - we must take into account all the forces on 
a subhalo, including the tidal force from nearby subhalos 
and the host halo core, in order to determine whether or 
not it truly is a bound system. Hence, previous algo- 
rithms, which do not do this, do not provide a complete 
picture. 

The final refinement to be made, therefore, is to ac- 
count for all the remaining forces on each particle in the 
subhalo. In order to do this we must now include the 
effect of the unbound particles on the potential energy 
of each subhalo during each iteration of the unbinding, 
and then calculate the tidal force on each particle in the 
subhalo due to the rest of the particles in the halo. The 
former is a simple change to the unbinding algorithm. 
Computing the tidal force on each subhalo particle, how- 
ever, is a more complex and costly procedure, especially 
in higher resolution halos. Since gravitational forces are 
linear, we can determine the total force on a particle as 
the sum of three terms: 

1. forces due to the particles in the unit considered, 

2. forces due to particles within the unit considered, 
but not bound to that unit, 

3. forces due to particles outside of the unit consid- 
ered. 

In this paper we describe a fast algorithm we have de- 
veloped in order to approximate the external tidal forces 
on a substructure. For the first time, we present a 'coher- 
ent' definition of a subhalo that accounts for all the forces 
it is subjected to, both internal and external, and thus 
identify all particles that are not just instantaneously 
bound to the substructure, but will remain so within 
the local dynamical timescale. We henceforth refer to 
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the method of identifying substructures used in previous 
studies as the standard subhalo definition. Essentially, 
these use the forces included in the first term above. We 
refer to the new approach outlined in this study as the 
coherent subhalo definition, allowing for the forces due 
to all three terms above. This is for the purpose of iden- 
tifying the groups of particles that remain together over 
many dynamical timescales, and thus correspond to the 
luminous galaxies they host. 

In Section [51 we summarize the main st eps of the al- 
gorith m common to both procedures (see IWeller et al.l 
(2005) for a complete description), and describe in detail 
the routine we have developed to calculate the tidal forces 
on each particle in a subhalo. We then demonstrate the 
impact of each stage of the standard and coherent sub- 
halo definitions on the substructure populations for a 
single test halo. In the last part of this section, we assess 
the accuracy of each method by tracking a small sample 
of subhalos over several timesteps, checking whether the 
particles identified as bound by each definition actually 
do stay with the subhalo over time. In Sec. [SJ we then 
compare the mass and radial distributions of substruc- 
tures for a large sample of halos analyzed using both 
definitions. We also look at the mass contained within 
subhalos of subhalos - or sub-substructures - in each 
sample to compare the subhalo populations of successive 
generations of subhalos in cluster-mass halos. 

2. DEFINITION OF SAMPLE 
2.1. The Identification of Halos and Subhalos 

In our preceding study (|Shaw et al.l I2006T ) . we pre- 
sented a comprehensive analysis of the physical char- 
acteristics of a large sample of cluster mass halos. 
In order to compare our results with previous work, 
we ado pted the standard s ubhalo definition, as out- 
lined in IWeller et al.l (|2005l ). This algorithm was ap- 
plied to a ACDM N-body simulation of 1024 3 parti- 
cles with box size 320/i~ 1 Mpc, particle mass (m p ) of 
2.54 x 10 9 /i _1 Mq and a spline kernel force softening 
length of e = 3.2/i _1 kpc. The simulation was evolved to 
z=0.05 using the Tree- Particle-Mesh (TPM) algorithm 
(|Bode fc Ostrikerl 12003). The cosmological parameters 
used include Sl m = 0.3, A = 0.7, and cts=0.95; outputs 
from this run have previously b een used to make predic- 
tions concerning strong lensing ( Wambs ganss "eTal][200l 
iHennawi et "all 120051: iDas & Ostrikeril2006f K We wish to 
ensure that all the halos are well resolved and that the 
overmerging problem is not in evidence. Thus we dis- 
card all halos with a virial mass less than 10,000 particles 
(« 3 x lO 13 ^-^/©). 

The method described in IWeller et all (|2005l ) starts 
with the hierarchical identification of structures. First, 
cluster mass halos are identified in the simulation box 
using a Friends-of-Friends routine with a linking length 
of b = 0.2n -1 / 3 , where n -1 / 3 is th e mean inter-particle 
separa tion. The Denmax routine of lBertschinger fc Gelbl 
(1991) is then run once on each FOF halo, using a high 
resolution smoothing length of 5e in order to identify the 
substructures. A family tree is then constructed by hi- 
erarchically associating the smallest mass subhalos with 
their lowest mass 'ancestor', so that each subhalo has 
only a single immediate parent. Those that consist of 
less than 30 particles are dissolved into their immediate 



parent. 

Up until this point the analysis is purely geometrical. 
Next, we must refine our identification of subhalos by 
determining their dynamically bound constituents. This 
is achieved in two stages: as a first approximation we 
calculate the center of mass velocity of each subhalo 
and remove those p articles that are st atistical outliers 
(see Section 2.2.1 in IWeller et~aTll2005l ). This step effi- 
ciently removes the most unbound particles in each sub- 
halo, thus allowing a more accurate determination of the 
subhalo center of mass velocity. We then complete the 
calculation exactly by iteratively identifying those parti- 
cles with a total energy greater than zero (in the center 
of mass frame of the subhalo) and removing the most 
energetic. At this point the standard and coherent sub- 
halo definitions diverge. In the former, once a particle is 
identified as 'unbound' (Eb in d > 0), it is removed entirely 
from future iterations of the calculation, and therefore its 
contribution to the gravitational potential energy of the 
substructure is neglected. This is the procedure that is 
adopted in the unbinding step by many publicly avail- 
able codes. However, in the coherent halo definition, wc 
include the contribution of the unbound particles within 
the subhalo. This has the effect of increasing the over- 
all binding energy of a subhalo, reducing the number of 
unbound particles that are identified. In either approach 
we dissolve a subhalo into its immediate parent if at any 
stage its mass drops below 30 particles. 

In both procedures, we next check to see if pairs of the 
immediate daughters of a parent cluster halo are bound. 
If this is the case, they form a hyper 'structure] the less 
massive of the two becomes a daughter of the more mas- 
sive structure. If a hyperstructure is found, we then check 
to see whether nearby particles previously not associated 
with either of the subhalos are bound to it. For the stan- 
dard method of analyzing halos, we finally remove all 
subhalos and any particles that are not bound to the en- 
tire structure and the procedure is complete. However, 
in our coherent subhalo definition, we now allow for the 
tidal forces on the subhalo due to the external mass in 
the cluster; a full description of how this is implemented 
is presented in the following section. Once this is done, 
we finally remove unbound particles and subhalos from 
the cluster, as in t he standard analysi s. 

As described in IWeller et al.l (|2005) , both procedures 
are both stable and largely independent of arbitrary pa- 
rameter choices. Henceforth, we refer to the most mas- 
sive structure in each cluster as the 'mother' or 'host' 
halo. We will quote the mass of the entire cluster in 
terms of its virial mass, M V i r . For a ACDM cosmology, 
it is conventional to define the virial mass M v ; r and radius 
i?vir as M vir = |7ri? 3 ir A c (2)p c (z), where p c is the crit- 
ical density of the universe, and the m ean over-density 
A c = 178ft m (z)°- 45 (|Lahav et al.lll991l ). In order to cal- 
culate the virial mass of each mother halo, we start at 
its density maximum and proceed outward until wc reach 
the virial radius, within which the mean over-density is 
A c ; we include all substructures with centers of mass 
within the virial radius. Hence we define a particle as 
being part of the halo if it either belongs solely to the 
mother halo or is part of a bound subhalo that has its 
center of mass within the virial radius of the mother. A 
flowchart of all the main steps in the coherent subhalo 
definition can be viewed in Fig. [TJ 
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Identify cluster mass halos in simulation volume 
using Friend s-of-Friends 



Identification of substructures within each halo 
using DENMAX 



Construct substructure family trees 



X 



Iteratively remove velocity outliers 



T 



Iteratively remove energetically unbound particles 



Search for hyperstructures and check again for 
nearby bound particles 



Remove tidally unbound particles from each 
subhalo 



Remove particles and substructures not bound to 
the entire assembly and those outwith R v i r 



Fig. 1. — Flowchart summarizing the main steps in the identifica- 
tion of subhalos — according to our new coherent subhalo definition 
- in cluster mass halos extracted from dissipationless cosmological 
N-body simulations. 

2.2. Calculation of external forces: the tidal 
approximation 

For a subhalo well separa ted from its mother halo, t he 
tidal radius Rt is given by (Bi nnev fe Tremaineill987t ) 



Rt = r d 



/ M(R t ) 
\3M(r d ) 



1/3 



(1) 



where r d is the distance between the density peak of the 
mother structure and the density peak of the substruc- 
ture, M(r d ) is the mass of all the particles within this ra- 
dius and M(R t ) is the mass of all the particles within the 
tidal radius (measured from the center of the subhalo). 
This formula simply asks if a particle feels a tidal force 
pulling it away from the substructure larger than the 
force pulling it to the center of the substructure. How- 
ever, this equation does not apply for most halos under 
consideration here, since they are well within the mother 
halo. We note that lKim fc P ark (200(f) have recently de- 
veloped a routine that uses the tidal radius to demarcate 
subhalo regions. 

Thus, in order to decide if a particle in a daughter halo 
will be tidally removed, we must first calculate the tidal 
force as the difference between the force on the particle 
and the force on the density peak of the daughter: 

Ftid = F p — Fd , (2) 

where for F p and Fd we sum up the force of all particles 
from the mother structure, including all other daughters, 
outside a radius Rt ■ Particles within the tidal radius are 
not contributing to the tidal (i.e. external) force. In order 
to decide if a particle is tidally bound we calculate the 
energy 



Etid = 



F t id(<')-vrci(t') dt' « F tid (t)-v rc i(i)r , (3) 




Fig. 2. — Schematic description of the quantities for the calcu- 
lation of the tidal radius. The big sphere represents the mother 
halo, and the small the daughter. The dark arrow in the middle is 
the vector to the density peak of the daughter. The box represents 
the face of the cube of length 2Rt which we use to calculate the 
approximate tidal forces. The light arrows illustrate three of the 
positions used to calculate the forces at the faces of the cube. 

where v re i is the velocity of the particle relative to the 
center of mass of the daughter, and r is a characteristic 
time. As the characteristic time, we use one fifth of the 
smaller period of either the particle's orbit around the 
center of the daughter halo or the daughter halo's or- 
bit around the center of the mother; this is the interval 
over which the tidal force would be approximately con- 
stant. This requirement ensures that there is little scope 
for varying r as a free parameter. Eqn. [3] thus gives an 
approximation for the energy gain of a particle due to 
the tidal force in the interval before this force changes 
substantially. If a particle has 



Em + E to t > 



(4) 



we declare the particle tidally unbound and move it to 
the associated mother. However with this criterion a 
small number of particles far from the substructure might 
be considered bound even when the tidal force is large, 
should its velocity by chance be perpendicular to the 
tidal force, leading Eqn. [3] to give a zero tidal energy. We 
avoid this problem if we include an additional criterion: 



1 tid 



W\ 



8 GM (|<5r|) m p 
> 3 



(5) 



where 5r is the position of the particle relative to the 
central density peak of the daughter, and M (\Sr\) is the 
mass within the daughter up to this radius. This is a re- 
laxed version of the tidal radius criterion in Eqn. [TJ and 
corresponds to the particle being at twice the tidal ra- 
dius if the two halos are spatially separated. We declare 
particles to be tidally unbound if they fulfil cither Eqn. [J] 
or Eqn. [5] 

A problem with the prescribed scheme above is that 
in order to calculate the tidal force for each particle we 
need to perform a costly N 2 operation. However, since 
an approximate correction for the tidal force will suffice, 
we calculate F t id with a linear approximation in the fol- 
lowing way. We consider a cube with its center at the 
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Fig. 3. — Schematic representation of the quantities required to 
establish if a daughter is tidally bound. 

density peak of the daughter and a side length of 2i? t as 
in Fig. [21 We calculate the tidal force exactly at the cen- 
ter of each of the six faces of the cube. We then perform 
a linear fit to calculate the tidal force at the center of the 
each of the faces with 



' tid 



A(r F 



(6) 



where is the peak density position of the daughter (rel- 
ative to the peak density position of the mother halo), 
and the entries of the 3x3 matrix A are obtained by a 
least squares fit of the exact tidal force at the center of the 
6 faces. This approximation corresponds to a divergence 
free field and hence to a quadrupole approximation, and 
works well when the tidal radius of the daughter is suffi- 
ciently small compared to the distance between the den- 
sity peak of the daughter and the mother. We use the 
approximation if i? t / |rd| < 0.5; otherwise we switch to 
an exact iV 2 calculation. We found that the median er- 
ror due to the linear tidal force approximation was about 
20%. 

The next step is to remove tidally unbound daughters 
from structures. From the point of view taken here, none 
of the daughters can be tidally removed from the entire 
cluster, because we are unable to calculate the tidal forces 
from the mass distribution surrounding it. So this step 
will not change the total amount of substructure in the 
halo, but it might reduce the size of hyperstructures. The 
criterion to determine if a daughter is tidally bound to a 
another substructure is 



Ar < R 



3M 



1/3 



(7) 



with M the mass of the mother of the substructure, tom 
the mass of the substructure, rrid the mass of the daugh- 
ter of the substructure, R the distance from the mother 
to the density peak of the substructure, and Ar the dis- 
tance between the substructure and the density peak of 
the daughter. These quantities are schematically repre- 
sented in Figure [3) 

2.3. Application to a Single Halo 



In order to compare the standard and coherent halo 
definitions, in Table [T] we give for each step the number 
of particles (N p ), the number of subhalos and the frac- 
tion of mass contained in substructure (/ s ) for one of the 
most massive halos in our sample, analyzed using both 
methods. The first line in the table gives the proper- 
ties of the halo at the point just prior to where the two 
methods diverge - after the initial geometrical identifi- 
cation of subhalos using DENMAX and havi ng removed 
the v elocity outliers from each subhalo (see IWeller et al.l 
2005). During the following two stages - the unbind- 
ing and identification of hyperstructures - over half the 
substructure is removed in both methods. However, at 
this point there is nearly 77% more substructure in the 
halo with the coherent analysis relative to the standard. 
This is because in the coherent analysis the unbound 
particles are retained in the potential calculation (but 
ignored in the standard analysis). Therefore, the gravi- 
tational potential experienced by each particle is greater, 
thus decreasing the number of particles identified as be- 
ing 'unbound'. However, once the tidally unbound par- 
ticles have been removed, the fraction of mass contained 
in substructure is only slightly greater for the coherent 
analysis. N p , f s and the number of subhalos all decrease 
further in the final stage as particles and subhalos not 
bound to the entire cluster, or those that are outwith 
the virial radius, are removed (see the last line of Table 
®. 

Fig. |4] displays projections of the density of substruc- 
ture particles in the halo after each step in Table [T] The 
greater mass of subhalos after the unbinding step in the 
coherent analysis relative to the standard analysis is ev- 
ident, as is the impact of the tidal step. 



2.4. Tracking Subhalo Mass 

We have described a new method for defining and mea- 
suring the bound mass of substructures in N-body simu- 
lations. However, in order to demonstrate decisively that 
our 'coherent' scheme is an improvement on the standard 
procedure, we must show that the particles it identifies 
as 'bound' are still associated with the subhalo at a later 
time, and that it is more accurate in doing so. Like- 
wise, we must demonstrate that the particles identified 
as 'unbound' are not associated with the subhalo at fu- 
ture timesteps. 

To perform this test, we use the output of a 256 3 
particle simulation (with particle mass m p — 3.15 x 
10 7 h' 1 M ) in a box of side-length 19.25 hr 1 Mpc 
and a (spline) softening length of 2.5h~ 1 kpc. We 
adopt a ACDM cosmology, with Qa = 0.733, H n = 
71 kms~ 1 Mpc~ 1 and erg = 0.84. The output of this sim- 
ulation was saved after every 10 particle-mesh timesteps 
(see Table[2|); henceforth, we refer to the time elapsed be- 
tween each saved output as being one 'timestep'. Each 
particle was allocated a unique identifier number at the 
beginning of the simulation, to enable us to track in- 
dividual particles over time. Although this simulation 
contained fewer particles that our main 1024 3 simulation, 
the smaller box size ensured that we could easily locate 
and follow a small number of well resolved cluster-mass 
halos. 

In each of the final six outputs of the simulation (from 
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Coherent Analysis 


Standard Analysis 


Step 


N p f s Nsubs 


N p f s Nsubs 


Dcnmax 

Subhalo unbinding and 
Hyperstructures 

Tidal Step 

Cluster unbinding and 
virial cut 


472659 34 8% 2185 
472659 8.4% 61 

472659 5.2% 58 
402411 4.5% 43 


472659 34 8% 2185 
472659 4.75% 56 

n/a 

402534 4.0 % 41 



TABLE 1 

Comparison between subhalo definitions of the substructure content of Halo 1 after the completion of each step in the 

respective procedures. 



z = 0.13 to 0, or a time period of 1.17 Gyrs) we iden- 
tified all structures of mass greater then 10 13 M Q using 
the procedure outlined in Sec. 12.11 From this sample 
we selected the most massive halo possessing a stable 
mass - one in which the mass fluctuates by less then 5% 
- over the time-period considered. The main physical 
properties of this halo can be viewed in Table El We also 
checked to make sure that the halo did not undergo any 
major mergers during this period. 

At each timestep, we applied both the standard and 
coherent subhalo definitions to identify the substructure 
populations of the halo. For each subhalo we identify five 
sets of particles. The set A contains all the particles geo- 
metrically associated with a subhalo, as identified by the 
DENMAX step (step 2 in Figure Q}, at time U (where 
the index i indicates the timestep) . The subset Ui con- 
tains all the particles that were associated with a subhalo 
by the DENMAX routine (i.e. that are in A) but were 
then removed during the unbinding stages. The subset 
Bi contains all particles that remained in the subhalo 
at the end of our entire routine - these are the 'bound' 
particles. Hence, for each subhalo we have five sets at 
each timestep; A, which is of course defined before the 
binding criteria are applied and so is the same for both 
the standard and coherent definitions, Ui and Bi as de- 
termined by the standard definition and Ui and Bi as 
determined by the coherent. The number of particles in 
any set Si is written N(S{). For example, the total num- 
ber of particles associated with a subhalo at a particular 
timestep is N(A) = N(Z7i) + N(Sj). 

To achieve our aim of following the bound constituents 
of a subhalo, we must be able to identify the same sub- 
halo at a later time. This is more complicated then sim- 
ply matching up particle id tags; over the course of time, 
subhalos may dissolve (by dropping below the 30 particle 
limit), merge, fragment, or temporarily move outwith the 
virial radius of the cluster or disappear into its core where 
DENMAX is unable to locate them. Hence, to achieve a 
one-to-one correspondence between subhalos at two dif- 
ferent times requires a clear definition of how to identify 
its descendant. To resolve this issue we developed a sim- 
ple matching algorithm that identifies the descendant of 
a subhalo when there might be some ambiguity. This 
routine operates by pairing the most massive progenitor 



with its most massive descendant, should there be more 
than one in either case. 

At each timestep considered, if, we look for the same 
subhalo at the subsequent timestep tj+j.. We do this for 
both the standard and coherent methods. If no descen- 
dant at the subsequent timestep is found for a particular 
subhalo for either method it is removed from the sample. 
This ensures a consistent comparison between the two 
methods (i.e. the same subhalos are being compared). 
Normally, a subhalo will only permanently disappear if 
it drops below the 30 particle minimum mass and is dis- 
solved into its immediate parent. For each subhalo, we 
are now able to calculate quantities such as N(AnA+i)> 
or the total number of particles that are found to be 
bound to the subhalo at both ti and £j+i. 

We now define three quantities (to be measured for 
each subhalo) that can be used to assess the accuracy of 
each criterion: 

• The probability P& that a particle in the subhalo 
is contained in the 'bound' set Bi but is not in 
the 'bound' set A+i at the subsequent timestep 
(written A+i): 



P, 



N(g, n (A+i)) 

N(A) 



(8) 



Note that B i+1 / Ui+i. This is because Ui + i 
includes only those particles that are geometri- 
cally associated to a subhalo at tj+i, but were 
removed from the subhalo during the unbinding 
stages. Bi + \ includes any particle in the entire halo 
that is not bound to the subhalo in question at U + i. 

• The probability P u that a particle in the subhalo 
is contained in the 'unbound' set Ui but is found in 
the 'bound' set A+i at the subsequent timestep: 



P„ 



N(c7 4 n Bi 



+1) 



N(A) 



(9) 



• The error on the bound mass determination for a 
subhalo, 



N(A)-N(A+iA) 

N(A) 



(10) 
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^lookback 

[Gyrs] 


M ha lo 
[fc _1 Af«,n] 


N sub (coh) 


N 3ub (std) 


35 


0.13 


1.17 


1.40 x 10 13 


59 


58 


36 


0.1 


0.92 


1.41 x 10 13 


56 


56 


37 


0.07 


0.66 


1.43 x 10 13 


56 


55 


38 


0.05 


0.48 


1.44 x 10 13 


57 


54 


39 


0.02 


0.19 


1.45 x 10 13 


60 


50 


40 





0.00 


1.46 x 10 13 


60 


52 



TABLE 2 

Mass and the number of subhalos (as identified by the standard (std) and coherent (coh) methods) in a halo followed 
over the final 6 steps of a 256 3 particle simulation (discussed in section 12.411 . 



or in words, the fraction of particles that were la- 
belled 'bound' at t t , minus the fraction of particles 
from the total set Di that actually were bound (i.e. 
are also found in Bi+i). 

The reason why we must write (Bi + \\Di) is that sub- 
halos sometimes 'sweep up' particles from the mother 
halo between timesteps. Of course, these particles can- 
not be accounted for by either criteria - we therefore 
ignore them in this analysis. 

If a criterion for identifying the bound mass compo- 
nent of subhalos in simulations claimed to operate per- 
fectly, we might expect to find = (particles that we 
say are bound, do actually stay bound), P u = (parti- 
cles that we say are unbound move away) and e m = 
(the number of particles identified as bound at ti is the 
same as the number identified as bound at tj+i, ignoring 
particles picked up along the way). In practice, within 
one timestep a subhalo may undergo a close encounter 
with another subhalo or the host halo core and there- 
fore lose particles due to physical processes that cannot 
be accounted for by applying dynamical criteria at each 
timestep. However, these effects will apply equally to 
subhalos identified by both methods, and thus will not 
affect our study of their relative accuracy. 

In Figures [5j[6] and [7] we plot the distribution of each of 
these quantities for the full sample of subhalos followed 
over a single timestep (i.e. i = 35-36, 36-37, etc). This is 
equivalent to 0.23 Gyrs, or 2.88 times the average char- 
acteristic dynamical time (r in Section [2~2]) . It is clear 
from the distributions of P u in Figure [5] that the coher- 
ent definition (dashed line) is less likely to mislabel as 
'bound' a particle that does not remain with a subhalo 
than the standard definition (solid line). The mean value 
of Pb for the standard sample (0.089 ± 0.006, where the 
error is the standard error in the mean) is 28% greater 
than that for the coherent sample (0.070 ±0.004). Hence, 
the standard routine is less successful at removing parti- 
cles that are not bound to the subhalo. However, both 
methods do equally well at correctly identifying the un- 
bound particles - the mean value of both distributions 
in Figure ©is (P u ) = 0.007. 

The distributions of e m in Figure [7J which are largely 



greater then zero, suggest that both methods overesti- 
mate the bound mass of a subhalo. However, it is clear 
that the error is significantly less for the coherent def- 
inition ((e m ) = 0.015 ± 0.002, dashed line/arrow) than 
the standard subhalo sample ((e m ) = 0.022 ±0.002, solid 
line/arrow). Again, this demonstrates that the standard 
routine is less accurate at identifying and removing un- 
bound particles and thus provides a less reliable measure 
of subhalo mass. 

We now investigate the accuracy of each method as a 
function of time (number of timesteps). This is achieved 
by also matching subhalos to their descendants between 
two and five timesteps later (equivalent to 0.473 — 1.17 
Gyrs, or 2.88 - 14.6 characteristic dynamical times) and 
determining the fraction of particles in each subhalo at 
ti that were correctly identified as being bound or un- 
bound. In Figures [5] and [H] we analyse e TO and Ph as a 
function of the time between outputs (At). In the upper 
panel of both plots, each point represents the mean val- 
ues (and the standard error in the mean) of e m and Ph 
respectively, for the sample of subhalos followed over 1-5 
timesteps. Hence, the left-most point in the upper panel 
of Figure [5] represents the mean value of Ph measured 
over a single timestep (as marked by the arrows in Fig- 
ure [5]), the second point over two timesteps, and so on. 
As usual, the dashed line represents the results obtained 
for the coherent sample; the solid line those obtained for 
the standard sample. 

As expected, Pb, the probability that a particle la- 
belled as 'bound' at time ti has left the subhalo at a 
later time ti + j (where j=l-5), increases as At increases 
for both methods (see Figure [5]) . After a single timestep 
(« 0.23 Gyrs), P b = 0.07 and 0.09 for the coherent and 
standard methods respectively. When measured over 
five timesteps (1.17 Gyrs), this increases to 0.25 and 
0.28. The lower panel, which charts the difference in the 
value of Pb obtained between the standard and coher- 
ent method, demonstrates that the coherent procedure 
becomes increasingly more accurate with respect to the 
standard with larger At. A similar result is demonstrated 
in Figure [U where we plot e(At). Gradually, particles 
identified as 'bound' by each method are removed from 
the subhalo by dynamical friction or tidal stripping, and 
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Fig. 5. — The distribution of P;, (see Section [Z4t for the coher- 
ent (dashed) and standard (solid) subhalo samples, measured over 
a single timestep. The arrows mark the mean values of each dis- 
tribution, respectively. It is clear that the coherent method is less 
likely to mislabel as 'bound' particles that infact leave the subhalo. 
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Fig. 4. — 2d projection of the substructure in one of the most 
massive halos in our sample after the completion of each stage of 
the analysis for both the standard and coherent halo definitions: 
after the denmax step and removal of the velocity outliers (a), the 
unbinding and identification of hyperstructures (b), the removal of 
tidally inbound particles (c) and, finally, having removed particles 
and subhalos unbound to the entire structure (d). 

so the error on the predicted subhalo mass at fcj grows 
with At. However, as discussed above, the increase in 
e and P& will be partly due to interactions with other 
subhalos or the cluster core within the time period over 
which they are evaluated and therefore is not solely due 
to inaccuracies in either method. Nevertheless, the er- 
rors for the standard method grow slightly more rapidly 
than for the coherent. 

Overall, although both methods do well at determin- 
ing those particles that are bound to a subhalo and those 
that will soon move away, we have clearly demonstrated 
that our new coherent method provides the more precise 
measure of subhalo mass. Whilst both the standard and 
coherent definitions overestimate subhalo masses, the lat- 
ter is less prone to mistakenly labelling 'unbound' parti- 
cles as 'bound'. It is important to note that the improve- 
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Fig. 6. — The distribution of P u (see Section [2.4| l for the coherent 
(dashed) and standard (solid) subhalo samples, measured over a 
single timestep. The arrow marks the mean values (which are 
the same) of both distributions. On average, the results indicate 
that, for both methods, less then 1% of the particles identified 
as 'unbound' are actually bound to (and thus remain with) the 
subhalo. 
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Fig. 7. — The distribution of e m - the erro r on the bound mass 
determination for each subhalo - (see Section l2.4H for the coherent 
(dashed) and standard (solid) subhalo samples, measured over a 
single timestep. The arrows mark the mean values of each distri- 
bution, respectively. The coherent method is clearly more accu- 
rate at measuring the bound mass of a subhalo than the standard 
method. 
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Fig. 8. — (upper) The mean value of Pi, obtained when measured 
over 1-5 timesteps (At = 0.23 - 1.17 Gyrs, or 2.88 - 14.6 char- 
acteristic dynamical times) for the standard (solid) and coherent 
(dashed) subhalo definitions. Error bars denote the standard error 
in the mean, (lower) The difference between the mean values of 
Pi, measured for the standard (std) and coherent (coh) methods as 
At increases. 




Fig. 9. — (upper) The mean value of e m obtained when measured 
over 1-5 timesteps (At = 0.23 - 1.17 Gyrs, or 2.88 - 14.6 charac- 
teristic dynamical times), for the standard (solid) and coherent 
(dashed) subhalo definitions. Error bars denote the standard error 
in the mean, (lower) The difference between the mean values of 
e m measured for the standard (std) and coherent (coh) methods 
as At increases. 

ment is evident despite the 20% error in the tidal forces 
due to the approximate scheme that we use to evaluate 
them (see Section |2~2|) . Additionally, at any instant, we 
are only able to determine those particles that will be 
removed by tidal forces within a single characteristic dy- 
namical timescale ((r) ~ 0.08 Gyrs = 0.35 timesteps), 
as described in Section |2~21 We therefore do not expect 
our coherent definition to be significantly more accurate 
then the standard scheme over timescales longer than a 
few dynamical times. 

Furthermore, to perform this analysis we picked a halo 
that contained a typical fraction of its mass in substruc- 
ture (~ 7%). If we had chosen a halo that contained 
several large substructures (i.e. a halo undergoing a ma- 
jor merger) then the tidal forces exerted on each subhalo 
would have been more substantial, requiring a larger cor- 
rection to their bound mass. 

3. COMPARISON BETWEEN METHODS FOR A LARGE 
HALO SAMPLE 

Both the standard and coherent subhalo definitions 
were applied to the 2000 most massive halos extracted 



from the large 1024 3 particle simulation desc ribed in Sec- 
tion 12.11 Using the procedure described in Sh aw et all 
(|200l we removed from each sample those halos that 
are not 'virialized', that is, have not yet reached a state 
of dynamical equilibrium; there then remains 1838 ha- 
los in each sample. In the following Section we compare 
the subhalo populations and their radial distributions for 
each sample, and finally we investigate the subhalo pop- 
ulations of subhalos. 

3.1. Substructure Populations of Halos 

Since the force and mass resolution of N-body simu- 
lations has reached a high enough level to enable the 
identification of subhalos within larger structures, several 
previous studies have attempted to measure the slope 
of the subhalo mass function (SMF). This is especially 
important if one wishes to quantify the relationship be- 
tween the properties of dark matter subhalos identified 
in simulati ons and observed properties of galaxies (see 
for exampl e Kravts ov et al.ll2~004atlTasitsiomi et al.ll200"4t 
IVale fc Os triker 2006) . Studies of small numbers of halos 
selected from low resolution simulations and resimulated 
at higher resolutions have suggested that dark matter 
halos are self-similar in terms of their subhalo popula- 
tions: low mass h alos appear like 'rescaled versions' of 
higher mass hal os (iMoore et al.lll999UGhigna et al.l l2000; 



iDe Lucia et~ld1l2004l : iReed et al.ll2005| ). However, recent 
studies of much larger samples have found evidence that 
this is not the case - higher mass halos contain more sub- 
structure, as they have formed more recently and have 
therefore had less time in which to disrupt their hosted 
subh alos (|Gao et al.ll2"00l iKang et al.ll2005t IShaw et all 
20061) . 

In Fig. [10] we compare the differential subhalo mass 
functions for each of our halos samples, quoting subhalo 
masses as a fraction of the cluster virial mass. The flat- 
tening off and eventual decline of each distribution for 
decreasing subhalo masses is due to the minimum sub- 
structure mass in our simulation (7.62 x lO 11 ^ -1 ./^©, 
or 30 particles). For subhalos with M su & > 0.01M V i r 
this effect appears not to be in evidence. Therefore, for 
subhalo masses greater than this, we fit a power-law - 
dN/dlog w (M) = (M sub /M vir ) a - to the distributions, 
obtaining a slope of —0.79 ± 0.04 for the coherent sam- 
ple (dashed lines) and —0.91 ±0.03 for the standard halo 
sample (solid lines). Typically, most published studies 
have found that a power-law fit to the slope of the sub- 
halo mass function resu l ts in values of a between -0.7 and 
1.1 (IMoore et al.lll999l: iGhigna et al.ll2000l: iHelmi et al 
Gao et all 12004 be Lucia et alj l2004USIiaw et al 
van den Bosch et al. 2005). Our results are there- 



2002; 



2006; 



fore well within the range measured by other authors. 

In general, there is good agreement between the two 
halo samples, although there is a greater number of very 
high mass subhalos (M su b > 0.2M V i r ) in the coherent 
sample, resulting in a shallower slope than measured 
for the standard sample. On average, these high mass 
subhalos tend to reside within slightly lower mass clus- 
ters ((M V ir) « 1.5 x 10 14 M Q ) than for the entire sam- 
ple ((M vir ) « 3 x 10 14 A/ Q ). They are also typically lo- 
cated near to the cluster core, where the density of back- 
ground particles - those bound only to the mother halo 
- is greater. In the coherent subhalo sample these ha- 
los will tend to be more massive than their counterparts 
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Fig. 10. — Mass distribution of all the subhalos in our coherent 
(dashed lines) and standard (solid) halo samples as a function of the 
ratio of subhalo to cluster halo mass. Straight lines are the power 
law fits to the slope outwards of the point indicated by the arrow, 
which is where we assume that the distribution is not affected by 
the minimum subhalo mass in our simulation (30 particles, or 7.62 X 
10 10 h- 1 M Q ). Slopes of the fit are -0.79 ± 0.04 and -0.91 ± 0.03 
for the coherent and standard distributions respectively. 




Fig. 11. — Maximum circular velocity distribution of all the 
subhalos in our coherent (dashed lines) and standard (solid) halo 
samples as a function of the ratio of subhalo to halo maximum 
circular velocity. Straight lines are the power law fits to the slope 
outwards of the point indicated by the arrow, which is where we 
assume that the distribution is not affected by the minimum sub- 
halo mass to which we can reliably measure subhalo maximum 
circular velocity (100 particles, or 2.54 X 10 Mq. Slopes of the fit 
are —3.13 ± 0.11 and —3.66 ± 0.30 for the coherent and standard 
distributions respectively. 

in the standard sample, as we include the contribution 
of the background particles to the binding energy of a 
subhalo. Furthermore, they are not subjected to strong 
tidal forces from the halo core because they are located 
in lower mass halos. 

Fig. QT] shows the circular velocity function (CVF) for 
subhalo maximum circular velocity as a fraction of halo 
maximum circular velocity, V smax /Vhmax- The maxi- 
mum circular velocity for both subhalos and halos was 
cal culated using the 3 parameter fitting formula proposed 
bv lStoehr! (pOOQ) (hereafter, STWS): 

log(-^) - -a{log(x)f (11) 

* max 

where x = r/r vmax and V max , r vmax and a are the max- 
imum circular velocity, the radius at which this is lo- 
cated, and the width of the profile, respectively. We 
have adopted this profile because the flexibility provided 



by the addition al free parameter com pared to the NFW 
density profile (|Navarro et al.l IT996f ). results in a more 
accurate measure of the maximum circular velocity of a 
halo. We find that the STWS formula provides a good 
fit to the circular velocity profile of subhalos of at least 
100 particles. Hence, in Fig. [TT]we only use subhalos of 
at least this mass. 

We fit a power-law to the distribution for V smax > 
OSVhmax, where it is less affected by the incompleteness 
effects of a minimum subhalo mass, obtaining a slope of 
—3.13 ±0.11 for the coherent halo sample (dashed lines) 
and — 3.66 ± 0.30 for the standard sample (solid lines). 
IReed et al.l £2005) measure a slope of -4 for their sam- 
ple of 16 high resolution subhalos, with significant scat- 
ter. The virial theorem gives the simple scaling relation, 
M oc V£ ax , where f3 = 3 for isolated halos. However, for 
subhalos having undergone mass loss, and because of a 
weak correlation between halo mass and concentration , 
tends to vary b e tween 3-4 (pWila- Reese et al.l Il999t 
iBullock et all 120011; lHavashi et all 120031 : IKravtsov et al l 
l2004bHShaw et aLll2006[ ). Hence, the relative values we 
obtain for the slopes of the SMF and CVF largely agree 
with what may be expected given the relationship be- 
tween M vir and V max . 

Overall, we again find a general good agreement be- 
tween the two samples, with evidence of a slightly higher 
number of high V sm ax subhalos in the coherent sample, 
resulting in the shallower slope. There is also a tiny 
fraction of subhalos with a greater maximum circular ve- 
locity than their host. These are very high mass sub- 
halos - M su i, w 0.38M vir - of almost equal mass to the 
mother halo itself. It is generally known that the maxi- 
mum circular velocity is a quantity that is fairly robust 
to changes in subha lo mass due to dynamic al friction and 
tidal stripping (e.g. IKravtsov et al.l (|2004bD ). Hence, it is 
not necessarily surprising that we should find high mass 
subhalos with V smax slightly greater than their hosts. 

In Fig. [TJ] we compare the distributions of the total 
fraction of mass contained in substructure (f s ) for both 
our halo samples. Aside from the slightly increased num- 
ber of very high mass substructures (relative to their host 
halo mass) discussed earlier, there is a very good agree- 
ment between the two definitions of substructure. The 
mean value of f s for the fully self-consistent definition of 
substructure halo sample is 0.089 ± 0.073, compared to 
0.082 ± 0.066 for the standard definition sample (errors 
correspond to one standard deviation) . These values also 
agree extremely well with those found in other studies 
(|De Lucia et al.l[200l iGao et al.ll200l iGill et al.ll2004bh 
who tend to find average values of f s in the range 8-10%. 

In Fig. [T3]we plot f s against halo mass for the coherent 
(top panel) and the standard (lower panel) halo samples. 
It is clear that for both samples that higher mass ha- 
los, on average, contain a higher fraction of their mass 
in substructure. To each distribution we fit a power law, 
obtaining similar slopes of 0.44 ±0.09 and 0.40 ±0.09 for 
the coherent and standard samples, respectively. This 
is an important result as it de monstrates that cluster- 
mass halos are not self-similar. IGao et all (|2004l ) found 
a similar result analyzing a r easonably large sampl e of ha - 
los over a wide mass range. Ivan den Bosch et~atl (2005) 
constructed a semi-analytical model to calculate the sub- 
halo mass function, adjusting the free parameters to 
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Fig. 12. — Distribution of the fraction of mass contained within 
substructure (f s ) for all the halos in our coherent (dashed lines) 
and standard (solid) samples. 
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Fig. 13. — Fraction of mass contained within substructure as a 
function of halo mass, for our coherent halo sample (upper panel) 
and our standard halo sample (lower). Solid lines are the best fit 
power law to each distribution, where we obtain slopes of 0.44±0.09 
and 0.40 ±0.09 for the coherent and standard samples respectively. 
We take M, as 8 X 10 14 M Q . 



match the s imulations oflGao et al.l (|2004f ) , iTormen et al.l 
(|2004f) and iDe Lucia et all (|2004D ; they found the slope 
and normalization of the mass function to be depen- 
dent on the ratio of halo mass to the charac teristic non- 
linea r mass scale. In our previous paper (I Shaw et al.l 
2006) we demonstrated that halos with a higher frac- 
tion of their mass contained in substructure typically 
also have a lower concentration, less spherical morphol- 
ogy and greater spin. We find the same results apply to 
the halos in our coherent sample. 

Overall, we find that accounting for all the forces on 
a substructure, both internal (including the contribution 
of the 'background' particles to the binding energy) and 
external (incorporating the tidally disruptive effect of the 
cluster mass), does not produce a significantly different 
subhalo population from the definition of substructure 
used in previous studies. In the following section, we an- 
alyze and discuss the radial distributions of substructures 
in our halo samples. 

3.2. Radial Distribution of Substructures 

The simplest assumption one may make is that the 
radial distribution of satellite galaxies in a cluster 



follows the distribution of the dark matter. However, 
a number of recent studies investigating the radial 
distribution of subhalos in dissipationless simulations 
of cluster mass halos have found that this distribu- 
tion is significantly less concentrated, and shallower 
in the inner r egion s, than that of th e dark matter 
(IGhigna et all 119981: I Colm et alj Il999t IGhigna et al l 
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2005t Faltcnbachcr & Dicmand 200$. iDiemand et~aTl 



(2004) performed convergence tests to show that the 
discrepancy between the dark matter density and sub- 
halo radial distributions is likely not due to numerical 
resolution issues in the inner regions of halos; rather, 
it is that the destructive forces to which a subhalo is 
subjected (dynamical friction, tidal shocks, merging) 
proceed more rapidly i n the inner regions than at the 
periphery of the system. INagai &: Kravtsovl (|2005l ) found 
in their simulations that whilst the radial distribution of 
subhalos was significantly less concentrated than that of 
the dark matter when selected using their present-day 
gravitationally bound mass, their radial distribution 
followed that of the dark matter much more closely 
when selected using their mass at the time of accretion. 
By performing simulations including gas dynamics, star 
formation and cooling, they also demonstrate that the 
radial distributions of subhalos selected by their stellar 
mass resulted in profiles similar in shape to that of 
the dark matter. As stars are located in the core of 
subhalos, and are therefore tightly bound, the stellar 
mass in a subhalo is a quantity that should not vary 
substantially over time. 

In order to compare the radial distributions of sub- 
structure in each of our samples, in the upper panel of 
Fig. [T3]we plot the spherically averaged fraction of mass 
in substructure within a given radius r, f s (< r/R V i r ). 
This is calculated by dividing the mass in substructures 
by the total mass inside spheres of successively increas- 
ing fractions of the virial radius; we then take the mean 
value of all the halos in each sample at each radii. The 
circles/dashed line show the coherent halo sample, and 
the diamonds/solid line show the standard sample. For 
reference, we also plot the mean fraction of the total 
halo mass, M{< r / R V i r ) / M (R V i r ) within each radius 
(crosses/dotted line). In Fig. [15] we plot the projected 
radial distribution of the mass contained in substructure; 
this is calculated by taking the mean of three orthogonal 
projections of each halo. In the upper panel we plot the 
substructure mass within a projected fraction of the halo 
virial radius, R/R V i r , as a fraction of halo virial mass, or 
M su b(< R)/M(R V i r ), for the coherent (dashed line) and 
standard (solid) halo samples. For reference, we also plot 
the total halo mass within R/R V i r as a fraction of M V i r 
(dotted line) . In the lower panel we plot the substructure 
mass within R/R V i r as a fraction of the total halo mass 
within R/R mr , or M sub (< R)/M DM {< R)- 

As found in previous studies, it is clear that the ra- 
dial distribution of substructure does not follow that 
of the dark matter in the inner regions of a cluster. 
We re-emphasize that the minimum subhalo mass in 
each sample is 30 particles or 7.62 x 10 10 /i _1 Mq and 
the minimum cluster virial mass is 10,000 particles or 
2.54 x 10 13 /i _1 Af o . Outward of 0.6i?„ ir , the cumula- 
tive fraction of mass in substructure for both samples re- 
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Fig. 14. — (upper panel) Cumulative fraction of mass in substruc- 
ture within spheres of radius r/R v i r , for our coherent (dashed line) 
and standard (solid line) halo samples. Dotted line denotes the to- 
tal mass contained with radius r as a fraction of the total halo 
mass at the virial radius, (lower) The fraction of mass in substruc- 
ture in radial bins, df/dr = (M au i,(r + dr) — M au ^(r)) / (Mtot(r + 
dr) — Mf„t(r)), for both halo samples. For reference we also plot 
(M to t(r + dr) - M to t(r))/M vir (dotted line). 

mains constant, indicating that in the outer regions the 
subhalos trace the overall dark matter distribution fairly 
well. Within this radius, the mass in subhalos relative 
to the dark matter decreases significantly. This is much 
more pronounced for the standard halo sample, which 
drops off much more steeply than the coherent sample. 
This is also shown in the lower panel of Fig. 1 141 where 
we plot the mean value of f s for each sample in radial 
bins. This is calculated by dividing the mass in sub- 
structure between r and r + dr by the total halo mass 
within the same spherical annulus. For reference we also 
plot (M to t{r + dr) — M tot (r)) / 'M vir at each radii (dotted 
line). Outwards of Q.5R V i r the results for each sample are 
very similar, but inside this radius there is significantly 
more substructure in the coherent halo sample. Thus in 
the inner regions of a halo, the binding effect of includ- 
ing the background (or unbound particles) located within 
a subhalo is on average slightly stronger than than the 
tidal forces from the cluster core and other surrounding 
substructures. 

3.3. Subhalos of subhalos 

When a halo is accreted onto a more massive host, it 
may bring with it its own internal hierarchy of substruc- 
tures. This is observed bv lGill et al.l (|2004bf ). who follow 
the temporal evolution of satellites in their simulations 
as they are accreted by their host. Assuming this inter- 
nal structure survives the initial impact of the disruptive 
environment within a cluster, we might expect that more 
massive halo s will contain severa l gen erations of sub- 
halos Both IZentner et ail (|2005l) and iTavlor fc Babul 

(2004) try to account for subhalos-of-subhalos during 
their construction of semi-analytic models for the sub- 
halo populations of dark matter halos. 

We have performed a preliminary investigation of the 
amount of mass that is contained in subhalo s that are not 
immed iate daughters of the mother halo. IWeller et al.l 

(2005) describe in detail how the family-tree style hier- 
archy of substructures in each halo is calculated, depend- 
ing on whether each subhalo is bound to the next most 
massive structure in the vicinity. Due to the resolution 
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Fig. 15. — (upper panel) Cumulative projected subhalo mass for 
the coherent (dashed lines) and standard (solid lines) halo samples 
as a fraction of halo virial mass, M su i,(< R)/M(R v i r ). Dotted line 
denotes the total halo mass contained within successively increas- 
ing fractions of the virial radius as a fraction of halo virial mass, 
Mtot(< R) I 'M (R v ir) . (lower panel) Cumulative projected mass in 
substructure in circles of radius R/R v i r as a fraction of the total 
halo mass within the circular region, M su ^(< r)/Mtot(< r). 

limits of our simulation, we are only able to resolve sub- 
structures down to a mass of approximately 10~ 3 — 10~ 4 
of the cluster halo virial mass. However, this is good 
enough to ensure that many of the cluster halos in each 
of our samples contain at least two generations of sub- 
structures. Indeed, several halos in our sample contain a 
small number of third generation subhalos. 

We have previously defined f s as the fraction of halo 
mass contained in substructure. We now also define f ss 
as the fraction of mass in each subhalo contained in sec- 
ond generation subhalos. Of course, the minimum sub- 
halo mass in our simulation, 7.62 x 10 10 /i _1 Mq or 30 
particles, inhibits our ability to identify sub-subhalos in 
subhalos more than it does to identify subhalos in ha- 
los. Therefore we only include in our sample subhalos 
with masses of at least 1% that of their host halo, or 
Msub > 0.01 M V i r , when calculating f s , and second gen- 
eration subhalos with M ss > 0.01 M su b when calculating 
/ ss . Additionally, we only analyze the sub-subhalo pop- 
ulations of subhalos that satisfy the first criterion, i.e. 
M su b > 0.01 M V i r . All halos and subhalos that contain 
no subhalos or sub-subhalos, respectively, are discarded. 
We therefore now have a consistent means of comparing 
the populations of different generations of subhalos. 

In the upper panel of Fig. [16]we compare the distribu- 
tion of f s for all the halos in our coherent sample (solid 
line) to the distribution of / ss for all the subhalos in the 
entire sample (dashed line) in logarithmic bins. In the 
lower panel, we plot the equivalent distributions for the 
standard halo sample. We normalize each distribution so 
that it is expressed as a fraction of the total number of 
halos with f s > for the former, and the number of ha- 
los with fss > for the latter. Hence, the distributions 
represent the probability of a halo containing f s of its 
mass in substructure, and the probability of a subhalo 
containing / ss of its mass in sub-substructure, assuming 
that both f s and / ss are greater than zero. 

The similarity between the distributions of f s and f ss 
- especially for the coherent sample - is striking. It 
suggests a self-similarity within the hierarchical struc- 
tures of cluster mass halos, with the 2nd generation of 
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Fig. 16. — Distribution of the fraction of the halo mass contained 
within substructure, f s (solid lines) and the fraction of subhalo 
mass contained within second generation subhalos, f sa (dashed 
lines) for our coherent (upper panel) and standard (lower panel) 
halo samples. Distributions are normalized by the total number of 
halos for which f s - and subhalos for which f ss — is greater than 
zero. 

subhalos distributed in the 1st generation, as the lat- 
ter is in the cluster as a whole. This result is both 
new and extremely interesting, especially considering 
that it has been demonstrated that halos are not self- 
similar - younger, higher mass halos tend to contain 
more substructure than their older, lower mass counter- 
parts (jGao et al.ll2004t [Piemand et aD2004t IShaw et all 
|2006| ). However, it may not entirely be surprising. Just 
as recently formed halos contain a higher fraction of their 
mass in substructure, recently accreted subhalos might 
be expected to be the same - they have only recently 
become exposed to the disruptive environment within a 
halo, and so still retain a significant number of their own 
subhalos. Older halos have had more time to disrupt the 
substructures they host, so older subhalos will have lost 
a higher fraction of their own sub-subhalos due to the 
same processes. We hope to investigate this further with 
higher resolution simulations, enabling the identification 
of several generations of subhalos, in the near future. 

4. DISCUSSION AND CONCLUSION 

In this paper we have presented a improved definition 
of subhalos in dissipationlcss dark matter N-body simu- 
lations, based on the coherent identification of their dy- 
namically bound constituents. Whereas previous meth- 
ods determine the energetically bound components of a 
subhalo while ignoring the contribution of particles in 
the halo that are not geometrically or dynamically as- 
sociated with it, our method accounts for all the forces, 
both internal and external, exerted on a subhalo. This 
is realized in two stages: first, when calculating the en- 
ergetically bound components of a subhalo, we include 
the contribution to the gravitational potential of those 
particles that were geometrically assigned to it, but that 
have been identified as unbound in previous iterations 
of the calculation. In previous studies, the contribution 
of unbound particles has been ignored. The effect of in- 
cluding their forces is to increase the binding energy, and 
therefore the mass of a subhalo. Second, we allow for the 
external forces exerted on each subhalo by the rest of the 
particles in the halo. This is achieved by approximating 
the tidal forces on each face of a cube surrounding each 



subhalo. We then calculate whether the tidal forces on 
each particle are sufficient to detach it from the subhalo 
within a characteristic time, discarding those for which 
this is the case. The effect of this stage therefore is to 
reduce the mass of subhalos - especially those near to 
the halo core or other dense conglomerates of particles. 

To demonstrate that our new coherent method of 
identifying subhalos results in a more accurate measure 
of their bound mass than the standard procedure, we 
tracked two samples of approximately 50 subhalos over 
six consecutive timesteps (from z = 0.13 to 0, or 1.17 
Gyrs), identified using the coherent and standard meth- 
ods respectively. By checking whether particles identified 
as bound to a subhalo (by either criterion) remained so, 
and whether nearby particles that were found not to be 
bound did indeed move away, we were able to assess how 
accurately each method identifies the bound mass of a 
subhalo. We demonstrated that the standard method is 
less successful at removing particles that are not bound 
to a subhalo. The average fraction of a subhalos mass 
that was identified as bound by the standard method but 
was found to have left the subhalo by the next timestep 
was, on average, a factor of 1.29 greater than that mea- 
sured for the coherent method (9% compared to 7%). 
As expected, the accuracy of both methods decreases 
when we compare subhalos to their descendants several 
timesteps ahead; after 5 timesteps (1.17 Gyrs or « 14 
characteristic dynamical times) 28% and 25% of the par- 
ticles identified as 'bound' by the standard and coherent 
methods respectively, were found to have left the sub- 
halo. 

We then applied both the standard method of iden- 
tifying subhalos and our new coherent definition to a 
sample of 1838 virialized halos extracted from a high res- 
olution cosmological simulation. We found the subhalo 
mass and maximum circular velocity distributions (rela- 
tive to the mass and maximum circular velocity of their 
hosts) for each sample to be very similar. The slope of 
the high mass end of the subhalo mass distribution was 
-0.79 ± 0.04 and -0.91 ± 0.03 for the coherent and stan- 
dard halo sample respectively. The slightly flatter slope 
for the coherent sample was due to a small increase in 
the number of very high mass subhalos, relative to their 
host. We obtained equivalent slopes of — 3.13± 0.11 and 
—3.66 ± 0.30 for the subhalo maximum circular veloc- 
ity distributions. The mean fraction of mass contained 
in substructure, / s , was 0.089 ± 0.074 for the coherent 
halo sample and 0.082 ± 0.066 for the standard sample. 
These results are all within the range of values found by 
previous authors. 

In agreement with iGao et all (|2004h , we find that f s 
tends to increase with halo mass for both halo samples. 
High mass halos (M„ r > M„) have formed more re- 
cently and have therefore have had less time to disrupt 
the su bhalos the halos t hey accreted. In our previous 
study (|Shaw et al.ll2006h . we demonstrated that - for 
our standard halo sample - halos with a high f s are typi- 
cally also less spherical, have a lower concentration and a 
higher angular momentum than halos with a lower frac- 
tion of their mass in substructure. We have verified that 
the same is true for our coherent halo sample. These 
results break the assumpti on that is adopte d in some 
semi-analytic models (e.g. lOguri &: Ledl2004[ ) that dark 
matter halos are self-similar, with galaxy mass halos re- 
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sembling 'miniature' cluster mass halos. Recently, using 
a prescription for modelling the evolution of subhalos un- 
der t he influence of dynam ical fraction and tida l strip- 
ping, IZentner et all (|2005f ) and iTavlor fc Babul (|2004f ) 
have developed semi-analytic models that correctly pre- 
dict the deviations from a self-similar scaling of subhalo 
abundance with halo mass. 

We then proceeded to investigate the radial distribu- 
tions of substructure in each sample. We found that 
within 60% of the halo virial radius the distribution of 
subhalos is significantly less concentrated than that of 
the dark matter particles. However, this effect was less 
pronounced for subhalos in the coherent halo sample, in- 
dicating that in the central regions, the binding effect of 
including the 'background particles' in the potential cal- 
culation is slightly stronger than the disruptive effect of 
the tidal forces. This increases the mass of the subhalos 
in the inner regions of a halo relative to their counter- 
parts in the standard halo sample. 

Overall, however, we found that the inclusion of the 
contribution of 'background' particles to the binding en- 
ergy of a subhalo was offset by the disruptive effect of the 
tidal forces from particles outside a subhalo. As these 
two effects, on average, tend to negate each other, we 
found that the subhalo populations were not substan- 
tially different to those found using the standard defini- 
tion of substructure. 

Finally, we performed a preliminary comparison of the 
subhalo population of halos to the second generation, 
or sub-subhalo, populations of subhalos. As part of our 
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subhalo identification scheme, we construct for each halo 
a 'family tree' of substructure, checking to see whether 
smaller subhalos are bound to more massive subhalos. 
Our results indicate that the distribution of 2nd genera- 
tion subhalos within their first generation hosts, mim- 
ics that of the distribution of subhalos in the overall 
halo. Hence, we find there to be a self-similar scaling 
between the populations of different generations of sub- 
halos within halos. This is a new and interesting result, 
especially given the lack of self-similarity in the overall 
subhalo populations of host halos of different mass. We 
intend to investigate this further in the near future, with 
higher resolution simulations, and thus with halos con- 
taining multiple generations of subhalos. 
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